During seismic survey of a subterranean region, seismic data are acquired typically by positioning a seismic source at a chosen shot location, and measuring the seismic reflections generated by the source using receivers placed at selected locations. The measured reflections are referred to as a single “shot record”. Many shot records are measured during a survey by moving the source and receivers to different locations and repeating the aforementioned process. The survey can then be used to perform Full-Wavefield Inversion, which uses the information contained in the shot records to determine physical properties of the subterranean region (e.g., speed of sound in the medium, density distribution, etc . . . ) Full-Wavefield Inversion is an iterative process, each iteration comprising the steps of forward modeling to create model data and objective function computation to measure the similarity between model and field data. Physical properties of the subsurface are adjusted at each iteration to ensure progressively better agreement between model and field data. Modification of subsurface properties needs to be carried out in such a way that known relationships between various properties are not violated. The update process typically generates several trial models, which may become unstable, leading to a “blow-up” (unbounded growth of the solution, until the numbers become so large that they can no longer be represented in a computer) of numerical simulations. Mathematically, a stable model corresponds to a positive semi-definite matrix of elastic constants (a matrix is positive semi-definite when all of its eigenvalues are non-negative), which enter as coefficients into the wave equation. The wave equation can be written in many different forms, depending on the level of physics that needs to be included in a simulation. For example, elastic propagation (a fairly general case) is described byρ(x)∂t2u(x,t)−∇·T(x,t)=g(x, t)T(x,t)=C(x):∇u(x,τ)≡cijkl(x)∂kul(x,τ)′where T is the stress tensor, x is a vector representing the three spatial coordinates, t is time, g is a source function, and cijpq, is a fourth-order tensor of elastic constants.
For convenience, cijpq is often mapped into a 6×6 matrix using Voight notation (Tsvankin (2005), see pg. 8):CIJ=cijkl, whereI=iδij+(1−δij)(9−i−j);J=kδkl+(1−δkl)(9−k−l);The CIJ matrix (or an equivalent matrix that represents the fourth-order tensor cijpq) needs to be positive definite (Helbig (1994), chapter 5).Current Technology
At a given iteration n of inversion, a model update usually involves computing a search direction sn (this is usually accomplished by computing the gradient of an objective function ƒ; often sn=−∇ƒ is used) and performing a “line search”, i.e., evaluating objective functions for various trial models which are created through a linear combination of a current model and the search direction:mn+1=mn+αsn 
The search direction is scaled by a “step size” α and added to the current model mn. The value of the scalar α that produces the best value of the objective function is selected and a new updated model is formed using this value. Sometimes (usually if the step size α is chosen to be too large) mn+1 may become physically infeasible and lead to a blow-up in numerical simulations. The blow-up may occur even if the model is unstable at only a few spatial locations. If this happens, one is forced to choose a different (usually smaller) step size, thus slowing down the inversion process.
Besides feasibility (stability) constraints described above, it may be appropriate to impose other constraints, e.g., require that all model parameters lie within a certain predetermined interval (“box constraints”). Such constraints are typically incorporated into the inversion process using penalty functions, Lagrange multipliers, or Projection Onto Convex Sets (POCS). The first two methods are appropriate when constraints are “soft”, i.e., can be violated at intermediate steps and need to be satisfied only at convergence. The last method, POCS, is appropriate for both soft and “hard” constraints (i.e., constraints which cannot be violated and need to be satisfied for all intermediate models). A conventional way of applying POCS to enforce soft constraints is to perform a projection at the end of the line search:mn+1=P[mn+αsn],where P is a projection operator. Hard constraints can be enforced in a similar manner:mn+1=mn+β(P[mn+asn]−mn).Fixing α, applying the projection operator P before the line search starts, and then performing a line search with 0<β<1 guarantees that all intermediate models will satisfy the desired constraint.